!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine hlld (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                           &
   ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, gm)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx
   real(8), intent(in) :: gm
   real(8), dimension(nx), intent(in) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l
   real(8), dimension(nx), intent(in) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r
   real(8), dimension(nx), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   integer(4) :: i, ip1, sgn
   real(8) :: s_m, s_l, s_r, s_1l, s_1r
   real(8) :: ro_1l, bx_1l, by_1l, bz_1l, en_1l
   real(8) :: ro_1r, bx_1r, by_1r, bz_1r, en_1r
   real(8) :: ro_2l, bx_2l, by_2l, bz_2l, en_2l
   real(8) :: ro_2r, bx_2r, by_2r, bz_2r, en_2r
   real(8) :: vx_l, vy_l, vz_l, pt_l, pr_l, v2_l, b2_l, s2_l, c2_l, cax2_l, cfx_l
   real(8) :: vx_r, vy_r, vz_r, pt_r, pr_r, v2_r, b2_r, s2_r, c2_r, cax2_r, cfx_r
   real(8) :: vx_1l, vy_1l, vz_1l, pt_1l
   real(8) :: vx_1r, vy_1r, vz_1r, pt_1r
   real(8) :: vx_2l, vy_2l, vz_2l, pt_2l
   real(8) :: vx_2r, vy_2r, vz_2r, pt_2r
   real(8) :: gmm1, ex1, ex2, ex3

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0

   do i = 1, nx - 1
!----------------------------------------------------------------------------------------------------------------------------------|
!     left vx, vy, vz, s, pt
!----------------------------------------------------------------------------------------------------------------------------------|
      vx_l = mx_l(i) / ro_l(i)
      vy_l = my_l(i) / ro_l(i)
      vz_l = mz_l(i) / ro_l(i)
      b2_l = bx_l(i) * bx_l(i) + by_l(i) * by_l(i) + bz_l(i) * bz_l(i)
      v2_l = vx_l * vx_l + vy_l * vy_l + vz_l * vz_l
      pr_l = (en_l(i) - v2_l * ro_l(i) / 2.0d0 - b2_l / 2.0d0) * gmm1
      c2_l = gm * pr_l / ro_l(i)
      s2_l = c2_l + b2_l / ro_l(i)
      cax2_l = bx_l(i) * bx_l(i) / ro_l(i)
      cfx_l = sqrt((s2_l + sqrt(s2_l * s2_l - 4.0d0 * c2_l * cax2_l)) / 2.0d0)

!----------------------------------------------------------------------------------------------------------------------------------|
!     right vx, vy, vz, s, pt
!----------------------------------------------------------------------------------------------------------------------------------|
      ip1 = i + 1
      vx_r = mx_r(ip1) / ro_r(ip1)
      vy_r = my_r(ip1) / ro_r(ip1)
      vz_r = mz_r(ip1) / ro_r(ip1)
      b2_r = bx_r(ip1) * bx_r(ip1) + by_r(ip1) * by_r(ip1) + bz_r(ip1) * bz_r(ip1)
      v2_r = vx_r * vx_r + vy_r * vy_r + vz_r * vz_r
      pr_r = (en_r(ip1) - v2_r * ro_r(ip1) / 2.0d0 - b2_r / 2.0d0) * gmm1
      c2_r = gm * pr_r / ro_r(ip1)
      s2_r = c2_r + b2_r / ro_r(ip1)
      cax2_r = bx_r(ip1) * bx_r(ip1) / ro_r(ip1)
      cfx_r = sqrt((s2_r + sqrt(s2_r * s2_r - 4.0d0 * c2_r * cax2_r)) / 2.0d0)

!----------------------------------------------------------------------------------------------------------------------------------|
!     left and right s, pt, sm
!----------------------------------------------------------------------------------------------------------------------------------|
      s_l = min (vx_l, vx_r) - max (cfx_l, cfx_r)
      s_r = max (vx_l, vx_r) + max (cfx_l, cfx_r)
      pt_l = pr_l + b2_l / 2.0d0
      pt_r = pr_r + b2_r / 2.0d0
      ex1 = (s_r - vx_r) * ro_r(ip1) - (s_l - vx_l) * ro_l(i) + 1d-12
      s_m = ((s_r - vx_r) * ro_r(ip1) * vx_r - (s_l - vx_l) * ro_l(i) * vx_l - pt_r + pt_l) / ex1
      pt_1l = ((s_r - vx_r) * ro_r(ip1) * pt_l - (s_l - vx_l) * ro_l(i) * pt_r +                                                     &
         ro_l(i) * ro_r(ip1) * (s_r - vx_r) * (s_l - vx_l) * (vx_r - vx_l)) / ex1
      pt_1r = pt_1l
      pt_2l = pt_1l
      pt_2r = pt_1l

!----------------------------------------------------------------------------------------------------------------------------------|
!     left var*
!----------------------------------------------------------------------------------------------------------------------------------|
      ro_1l = ro_l(i) * (s_l - vx_l) / (s_l - s_m + 1d-12)
      vx_1l = s_m
      ex1 = ro_l(i) * (s_l - vx_l) * (s_l - s_m) - bx_l(i) * bx_l(i) + 1d-12
      ex2 = (s_m - vx_l) / ex1
      vy_1l = vy_l - bx_l(i) * by_l(i) * ex2
      vz_1l = vz_l - bx_l(i) * bz_l(i) * ex2
      bx_1l = bx_l(i)
      ex2 = (ro_l(i) * (s_l - vx_l) ** 2 - bx_l(i) * bx_l(i)) / ex1
      by_1l = by_l(i) * ex2
      bz_1l = bz_l(i) * ex2
      ex1 = vx_l * bx_l(i) + vy_l * by_l(i) + vz_l * bz_l(i) - vx_1l * bx_1l - vy_1l * by_1l - vz_1l * bz_1l
      en_1l = ((s_l - vx_l) * en_l(i) - pt_l * vx_l + pt_1l * s_m + bx_l(i) * ex1) / (s_l - s_m + 1d-12)

!----------------------------------------------------------------------------------------------------------------------------------|
!     right var*
!----------------------------------------------------------------------------------------------------------------------------------|
      ro_1r = ro_r(ip1) * (s_r - vx_r) / (s_r - s_m + 1d-12)
      vx_1r = s_m
      ex1 = ro_r(ip1) * (s_r - vx_r) * (s_r - s_m) - bx_r(ip1) * bx_r(ip1) + 1d-12
      ex2 = (s_m - vx_r) / ex1
      vy_1r = vy_r - bx_r(ip1) * by_r(ip1) * ex2
      vz_1r = vz_r - bx_r(ip1) * bz_r(ip1) * ex2
      bx_1r = bx_r(ip1)
      ex2 = (ro_r(ip1) * (s_r - vx_r) ** 2 - bx_r(ip1) * bx_r(ip1)) / ex1
      by_1r = by_r(ip1) * ex2
      bz_1r = bz_r(ip1) * ex2
      ex1 = vx_r * bx_r(ip1) + vy_r * by_r(ip1) + vz_r * bz_r(ip1) - vx_1r * bx_1r - vy_1r * by_1r - vz_1r * bz_1r
      en_1r = ((s_r - vx_r) * en_r(ip1) - pt_r * vx_r + pt_1r * s_m + bx_r(ip1) * ex1) / (s_r - s_m + 1d-12)

!----------------------------------------------------------------------------------------------------------------------------------|
!     left and right s*
!----------------------------------------------------------------------------------------------------------------------------------|
      s_1l = s_m - abs (bx_l(i)) / sqrt (ro_1l)
      s_1r = s_m + abs (bx_r(ip1)) / sqrt (ro_1r)

!----------------------------------------------------------------------------------------------------------------------------------|
!     left and right var**
!----------------------------------------------------------------------------------------------------------------------------------|
      if (bx_l(i) .ge. 0) then
         sgn = 1
      else
         sgn = -1
      endif

      ex1 = sqrt (ro_1l)
      ex2 = sqrt (ro_1r)
      ex3 = ex1 + ex2
      ro_2l = ro_1l
      ro_2r = ro_1r
      vx_2l = s_m
      vx_2r = s_m
      vy_2l = (ex1 * vy_1l + ex2 * vy_1r + (by_1r - by_1l) * sgn) / ex3
      vy_2r = vy_2l
      vz_2l = (ex1 * vz_1l + ex2 * vz_1r + (bz_1r - bz_1l) * sgn) / ex3
      vz_2r = vz_2l
      bx_2l = bx_l(i)
      bx_2r = bx_r(ip1)
      by_2l = (ex1 * by_1r + ex2 * by_1l + ex1 * ex2 * (vy_1r - vy_1l) * sgn) / ex3
      by_2r = by_2l
      bz_2l = (ex1 * bz_1r + ex2 * bz_1l + ex1 * ex2 * (vz_1r - vz_1l) * sgn) / ex3
      bz_2r = bz_2l
      ex3 = vx_2l * bx_2l + vy_2l * by_2l + vz_2l * bz_2l
      en_2l = en_1l - ex1 * (vx_1l * bx_1l + vy_1l * by_1l + vz_1l * bz_1l - ex3) * sgn
      en_2r = en_1r + ex2 * (vx_1r * bx_1r + vy_1r * by_1r + vz_1r * bz_1r - ex3) * sgn

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux
!----------------------------------------------------------------------------------------------------------------------------------|
      if (s_l .gt. 0.0d0) then
         fro(i) = mx_l(i)
         fmx(i) = vx_l * mx_l(i) + pr_l + b2_l / 2.0d0 - bx_l(i) * bx_l(i)
         fmy(i) = vx_l * my_l(i) - bx_l(i) * by_l(i)
         fmz(i) = vx_l * mz_l(i) - bx_l(i) * bz_l(i)
         fbx(i) = 0.0d0
         fby(i) = vx_l * by_l(i) - bx_l(i) * vy_l
         fbz(i) = vx_l * bz_l(i) - bx_l(i) * vz_l
         fen(i) = (en_l(i) + pr_l + b2_l / 2.0d0) * vx_l - (bx_l(i) * vx_l + by_l(i) * vy_l + bz_l(i) * vz_l) * bx_l(i)
      elseif (s_l .le. 0.0d0 .and. s_1l .gt. 0.0d0) then
         fro(i) = mx_l(i)
         fmx(i) = vx_l * mx_l(i) + pr_l + b2_l / 2.0d0 - bx_l(i) * bx_l(i)
         fmy(i) = vx_l * my_l(i) - bx_l(i) * by_l(i)
         fmz(i) = vx_l * mz_l(i) - bx_l(i) * bz_l(i)
         fbx(i) = 0.0d0
         fby(i) = vx_l * by_l(i) - bx_l(i) * vy_l
         fbz(i) = vx_l * bz_l(i) - bx_l(i) * vz_l
         fen(i) = (en_l(i) + pr_l + b2_l / 2.0d0) * vx_l - (bx_l(i) * vx_l + by_l(i) * vy_l + bz_l(i) * vz_l) * bx_l(i)
         fro(i) = fro(i) + s_l * (ro_1l - ro_l(i))
         fmx(i) = fmx(i) + s_l * (vx_1l * ro_1l - mx_l(i))
         fmy(i) = fmy(i) + s_l * (vy_1l * ro_1l - my_l(i))
         fmz(i) = fmz(i) + s_l * (vz_1l * ro_1l - mz_l(i))
         fbx(i) = fbx(i) + s_l * (bx_1l - bx_l(i))
         fby(i) = fby(i) + s_l * (by_1l - by_l(i))
         fbz(i) = fbz(i) + s_l * (bz_1l - bz_l(i))
         fen(i) = fen(i) + s_l * (en_1l - en_l(i))
      elseif (s_1l .le. 0.0d0 .and. s_m .gt. 0.0d0) then
         fro(i) = mx_l(i)
         fmx(i) = vx_l * mx_l(i) + pr_l + b2_l / 2.0d0 - bx_l(i) * bx_l(i)
         fmy(i) = vx_l * my_l(i) - bx_l(i) * by_l(i)
         fmz(i) = vx_l * mz_l(i) - bx_l(i) * bz_l(i)
         fbx(i) = 0.0d0
         fby(i) = vx_l * by_l(i) - bx_l(i) * vy_l
         fbz(i) = vx_l * bz_l(i) - bx_l(i) * vz_l
         fen(i) = (en_l(i) + pr_l + b2_l / 2.0d0) * vx_l - (bx_l(i) * vx_l + by_l(i) * vy_l + bz_l(i) * vz_l) * bx_l(i)
         fro(i) = fro(i) + s_1l * (ro_2l - ro_1l) + s_l * (ro_1l - ro_l(i))
         fmx(i) = fmx(i) + s_1l * (vx_2l * ro_2l - vx_1l * ro_1l) + s_l * (vx_1l * ro_1l - mx_l(i))
         fmy(i) = fmy(i) + s_1l * (vy_2l * ro_2l - vy_1l * ro_1l) + s_l * (vy_1l * ro_1l - my_l(i))
         fmz(i) = fmz(i) + s_1l * (vz_2l * ro_2l - vz_1l * ro_1l) + s_l * (vz_1l * ro_1l - mz_l(i))
         fbx(i) = fbx(i) + s_1l * (bx_2l - bx_1l) + s_l * (bx_1l - bx_l(i))
         fby(i) = fby(i) + s_1l * (by_2l - by_1l) + s_l * (by_1l - by_l(i))
         fbz(i) = fbz(i) + s_1l * (bz_2l - bz_1l) + s_l * (bz_1l - bz_l(i))
         fen(i) = fen(i) + s_1l * (en_2l - en_1l) + s_l * (en_1l - en_l(i))
      elseif (s_1r .gt. 0.0d0 .and. s_m .le. 0.0d0) then
         fro(i) = mx_r(ip1)
         fmx(i) = vx_r * mx_r(ip1) + pr_r + b2_r / 2.0d0 - bx_r(ip1) * bx_r(ip1)
         fmy(i) = vx_r * my_r(ip1) - bx_r(ip1) * by_r(ip1)
         fmz(i) = vx_r * mz_r(ip1) - bx_r(ip1) * bz_r(ip1)
         fbx(i) = 0.0d0
         fby(i) = vx_r * by_r(ip1) - bx_r(ip1) * vy_r
         fbz(i) = vx_r * bz_r(ip1) - bx_r(ip1) * vz_r
         fen(i) = (en_r(ip1) + pr_r + b2_r / 2.0d0) * vx_r - (bx_r(ip1) * vx_r + by_r(ip1) * vy_r + bz_r(ip1) * vz_r) * bx_r(ip1)
         fro(i) = fro(i) + s_1r * (ro_2r - ro_1r) + s_r * (ro_1r - ro_r(ip1))
         fmx(i) = fmx(i) + s_1r * (vx_2r * ro_2r - vx_1r * ro_1r) + s_r * (vx_1r * ro_1r - mx_r(ip1))
         fmy(i) = fmy(i) + s_1r * (vy_2r * ro_2r - vy_1r * ro_1r) + s_r * (vy_1r * ro_1r - my_r(ip1))
         fmz(i) = fmz(i) + s_1r * (vz_2r * ro_2r - vz_1r * ro_1r) + s_r * (vz_1r * ro_1r - mz_r(ip1))
         fbx(i) = fbx(i) + s_1r * (bx_2r - bx_1r) + s_r * (bx_1r - bx_r(ip1))
         fby(i) = fby(i) + s_1r * (by_2r - by_1r) + s_r * (by_1r - by_r(ip1))
         fbz(i) = fbz(i) + s_1r * (bz_2r - bz_1r) + s_r * (bz_1r - bz_r(ip1))
         fen(i) = fen(i) + s_1r * (en_2r - en_1r) + s_r * (en_1r - en_r(ip1))
      elseif (s_r .gt. 0.0d0 .and. s_1r .le. 0.0d0) then
         fro(i) = mx_r(ip1)
         fmx(i) = vx_r * mx_r(ip1) + pr_r + b2_r / 2.0d0 - bx_r(ip1) * bx_r(ip1)
         fmy(i) = vx_r * my_r(ip1) - bx_r(ip1) * by_r(ip1)
         fmz(i) = vx_r * mz_r(ip1) - bx_r(ip1) * bz_r(ip1)
         fbx(i) = 0.0d0
         fby(i) = vx_r * by_r(ip1) - bx_r(ip1) * vy_r
         fbz(i) = vx_r * bz_r(ip1) - bx_r(ip1) * vz_r
         fen(i) = (en_r(ip1) + pr_r + b2_r / 2.0d0) * vx_r - (bx_r(ip1) * vx_r + by_r(ip1) * vy_r + bz_r(ip1) * vz_r) * bx_r(ip1)
         fro(i) = fro(i) + s_r * (ro_1r - ro_r(ip1))
         fmx(i) = fmx(i) + s_r * (vx_1r * ro_1r - mx_r(ip1))
         fmy(i) = fmy(i) + s_r * (vy_1r * ro_1r - my_r(ip1))
         fmz(i) = fmz(i) + s_r * (vz_1r * ro_1r - mz_r(ip1))
         fbx(i) = fbx(i) + s_r * (bx_1r - bx_r(ip1))
         fby(i) = fby(i) + s_r * (by_1r - by_r(ip1))
         fbz(i) = fbz(i) + s_r * (bz_1r - bz_r(ip1))
         fen(i) = fen(i) + s_r * (en_1r - en_r(ip1))
      elseif (s_r .le. 0.0d0) then
         fro(i) = mx_r(ip1)
         fmx(i) = vx_r * mx_r(ip1) + pr_r + b2_r / 2.0d0 - bx_r(ip1) * bx_r(ip1)
         fmy(i) = vx_r * my_r(ip1) - bx_r(ip1) * by_r(ip1)
         fmz(i) = vx_r * mz_r(ip1) - bx_r(ip1) * bz_r(ip1)
         fbx(i) = 0.0d0
         fby(i) = vx_r * by_r(ip1) - bx_r(ip1) * vy_r
         fbz(i) = vx_r * bz_r(ip1) - bx_r(ip1) * vz_r
         fen(i) = (en_r(ip1) + pr_r + b2_r / 2.0d0) * vx_r - (bx_r(ip1) * vx_r + by_r(ip1) * vy_r + bz_r(ip1) * vz_r) * bx_r(ip1)
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine hlld
